Effect of surface roughness on rate-dependent slip in simple fluids 
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Molecular dynamics simulations are used to investigate the influence of molecular- 
scale surface roughness on the slip behavior in thin liquid films. The slip length 
increases almost linearly with the shear rate for atomically smooth rigid walls and 
incommensurate structures of the liquid/solid interface. The thermal fluctuations of 
the wall atoms lead to an effective surface roughness, which makes the slip length 
weakly dependent on the shear rate. With increasing the elastic stiffness of the wall, 
the surface roughness smoothes out and the strong rate dependence is restored again. 
Both periodically and randomly corrugated rigid surfaces reduce the slip length and 
its shear rate dependence. 

PACS numbers: 68.08.-p, 83.50.Rp, 47.61.-k, 83.10.Rs 



I. INTRODUCTION 

The description of the fluid flow in conflned geometry requires specification of the bound- 
ary condition for the fiuid velocity at the solid wall. Usually the fiuid is assumed to be 
immobile at the boundary. Although this assumption is successful in describing fiuid fiow 
on macroscopic length scales, it needs a revision for the microscopic scales due to possible 
slip of the fiuid relatively to the wall I'l. The existence of liquid slip at the solid surfaces was 

nnn 

established in many experiments on the pressure driven fiow in narrow capillaries [2|, |3|, |J| 
and drainage of thin liquid films in the surface force apparatus . The most popular 

Navier model relates the fiuid slip velocity to the interfacial shear rate by introducing the 
slip length, which is assumed to be rate- independent. The slip length is defined as a distance 
from the boundary where the linearly extrapolated fiuid velocity profile vanishes. Typical 
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values of the slip length inferred from the experiments on fluids confined between smooth 
hydrophobic surfaces is of the order of ten nanometers j^, Q, y, y, Q. Despite the large 
amount of experimental data on the slip length HI], the underlying molecular mechanisms 
leading to slip are still poorly understood because it is very difficult to resolve the fluid 
velocity profile in the region near the liquid/solid interface at these length scales. 

Over the last twenty years, molecular dynamics (MD) simulations were extensively used 
to investigate the correlation between the structure o: 
ically smooth surfaces and slip boundary conditions 
advantage of the MD method is that the fiuid velocity profile can be resolved at the molecu- 
lar level and no assumptions about the slip velocity at the interface are required. The main 
factors affecting slip for atomically smooth surfaces are the wall-fluid interaction, the degree 
of commensurability of liquid and solid structures at the interface, and diffusion of fluid 
molecules near the wall. The slip length was found to correlate inversely with the wall-fluid 
interaction energy and the amount of structure induced in the flrst fluid layer by the periodic 
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surface potential p^. For weak wall- fluid interactions and smooth surfaces, the slip length is 



proportional to the collective relaxation coefficient of the fiuid molecules near the wall [20]. 
The thermal fiuctuations of the wall atoms under the strong harmonic potential reduce the 
degree of the in-plane fiuid ordering and result in larger values of the slip length 12|. On 
the other hand, an excessive penetration of the wall atoms into the fiuid phase reduces the 
slip velocity for soft thermal walls 21| . Nevertheless, the effect of thermal surface roughness 
on the slip length in the shear-rate- dependent regime was not systematically explored even 
for atomically smooth walls. 



In the original MD study by Thompson and Troian 22] on boundary driven shear flow 
of simple fluids past atomically smooth rigid walls, the slip length was found to increase 
nonlinearly with the shear rate for weak wall-fluid interactions. A similar dynamic behavior 
of the slip length has been reported in thin polymer fllms Q]. The rate-dependent slip was 
also observed for the planar Poiseuille flow of simple fluids conflned between hydrophobic 



surfaces with variable size of the wall atoms 
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25l ]. The variation of the slip length (from 



negative to positive values) with increasing shear rate for hydrophilic surfaces [24] can be 
well described by the power law function proposed in Ref. 22|. In the recent paper 26| . 
we have reported a gradual transition in the shear rate dependence of the slip length, from 
linear to highly nonlinear function with pronounced upward curvature, by decreasing the 
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strength of the wall-fluid interaction. Remarkably, in a wide range of shear rates and surface 
energies, the slip length is well fltted by a power law function of a single variable, which 
is a combination of the structure factor, contact density, and temperature of the flrst fluid 
layer. One of the goals of the present study is to investigate how the rate-dependent slip is 
affected by the presence of the molecular-scale surface roughness. 

Molecular scale simulations of simple 27|, |28|] and polymeric 29|, |30| fluids (as well as 



recent experiments 
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331] ) have shown that the slip length is reduced in the presence of 



the surface roughness. The effect is enhanced for smaller wavelengths and larger amplitudes 



of the surface corrugation [28|. 



At low shear rates, the reduction of the effective slip 



length is caused by the local curvature of the fluid flow above macroscopic surface corruga- 



tions [28|, 
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35( 1 or by more efficient trapping of the fluid molecules by atomic-scale surface 



inhomogeneities 
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36|. The analysis of more complex systems with combined 



effects of surface roughness and rate dependency poses certain difficulties in the interpre- 
tation of the experimental results because the exact dependence of the local slip length on 
shear rate is often not known. 

In this paper, we explore the influence of molecular-scale surface roughness on the slip 
behavior in a flow of simple fluids driven by a constant force. We will show that the 
functional form of the rate-dependent slip length is considerably modifled by the presence 
of the thermal, random and periodic surface roughness. The growth of the slip length 
with increasing shear rate, which is observed for atomically smooth rigid walls, is strongly 
reduced by periodic and random surface roughness. Soft thermal walls produce very weak 
rate dependence of the slip length, while the linear behavior is restored for stiffer walls. 

The paper is organized as follows. The details of molecular dynamics simulations are 
described in the next section. Results for the shear rate dependence of the slip length on the 
thermal and random surface roughness are presented in Section IIIII The effect of periodic 
wall roughness on the rate-dependent slip is discussed in Section IIVI The summary is given 
in the last section. 



II. MD SIMULATION MODEL 



The simulation setup consists of = 3456 fluid molecules conflned between two stationary 
atomistic walls parallel to the xy plane. The molecules interact through the truncated 
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Lennard- Jones (LJ) potential 

VUr)^,._[(^f-{^y], (1) 

where e and a represent the energy and length scales of the fluid phase with density 
p = 0.81 cr^'^. The interaction between wall atoms and fluid molecules is also modeled by 
the LJ potential with the energy e^f and length scale cr^f measured in units of e and a. 
In all our simulations, wall atoms do not interact with each other and o"„f = a. The cutoff 
distance is set to rc = 2.5(j for fiuid-fiuid and wall-fluid interactions. 

The steady-state flow was induced by a constant force fx in the x direction, which acted on 
each fluid molecule. The heat exchange with external reservoir was controlled by a Langevin 
thermostat with a random, uncorrelated force and a friction term, which is proportional to 
the velocity of the fluid molecule. The value of the friction coefficient r = 1.0r~^ is small 
enough not to affect significantly the self-diffusion coefficient of the fluid molecules 37 



m. 



The Lan gevi n thermostat was applied only along the y axis to avoid a bias in the flow 



direction 12] . The equations of motion for a fluid monomer of mass m are given by 

dxi 



-^^ = -Ei?'+f- (2) 



myi + mVyi = - + fi , (3) 

-^^ = -E§' (4) 

where fi is a randomly distributed force with (/j(t))=0 and variance 



TABLE I: Root mean-square displacement, {6u'^) = 3kBT/K, divided by the nearest-neighbor dis- 
tance d = 0.8a and the typical oscillation time of the wall atoms tethered about their equilibrium 
lattice positions as a function of the spring stiffness for = 4m. 
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{fi{0)fj(t)) = 2mkBTTS(t)Sij obtained from the fluctuation-dissipation theorem. The 
temperature of the thermostat is fixed to T = Lle/ks, where ks is the Boltzmann constant. 
The equations of motion of the fluid molecules and wall atoms are integrated using the 



39( 1 with a time step At = 0.002 r, where r= ^Jma'^/e is 



fifth-order gear-predictor method 
the LJ time. 

The wall atoms were allowed to oscillate about their equilibrium lattice sites under the 
harmonic potential Vsp = \ nr"^. The spring stiffness k was chosen so that the ratio of the 
root mean-square displacement of the wall atoms and their nearest-neighbor distance was 



less than the Lindemann criterion for melting, a/ {6u'^)/d < 0.15, e.g. see Ref. 40|. At the 
same time, the parameter k should be small enough so that the dynamics of the wall atoms 
can be accurately resolved with the MD integration time step. The mass of the wall atoms 
rriw was chosen to be four times that of the fluid molecule to make their oscillation times 
comparable. The Langevin thermostat was also applied to all three components of the wall 
atoms equations of motion, e.g. for the ct component 

Xi + m^Txi = - ^ - + fi , (5) 

where the sum is taken over the fluid molecules within the cutoff radius Vc = 2.5 cr. Mean 
displacement of the wall atoms and their typical oscillation time are summarized in Table [T] 
for the values of the spring constant considered in this study. 

Each wall is constructed out of 648 atoms distributed between two (111) planes of the 
face-centered cubic (fee) lattice with density pw = 2.73 a~^. The walls are separated by a 
distance 24.58 a along the z axis. The lateral dimensions of the computational domain in 
the xy plane are measured as Lj. = 25.03(t and Ly = 7. 22 a. Periodic boundary conditions 
are imposed along the x and y directions. After an equilibration period of at least 2 x lO^r, 
the fluid velocity and density profiles were computed for a time interval up to 2 x lO^r 
within bins of thickness Az = 0.2a and A2; = 0.01cr, respectively 26| . The shear viscosity, 
fi = (2.2 ± 0.2) £Ta~^, remained independent of the external force |23l. |26|. 
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FIG. 1: (Color online) Averaged fluid density profiles near the thermal K = 600e/(T^ (o) and rigid 
(o) walls for e„f/e = 0.9 and fx = 0.001 e/o". Left vertical axis coincides with the location of the 
liquid/solid interface at z = — 11.79 o". 

III. RESULTS FOR THE THERMAL WALLS 
A. Fluid density and velocity profiles 



Examples of the averaged fluid density proflles near the thermal and rigid walls are 
presented in Fig. [T] for a small value of the external force fx = 0.001 e/cr. The flrst peak in 
the density proflle is slightly broader for the thermal walls because the fluid molecules can 
move closer to the fee lattice plane due to flnite spring stiffness of the wall atoms. The 
maximum value of the flrst peak deflnes a contact density, which is larger for the rigid walls 
because of the higher in-plane fluid ordering. In both cases, the fluid density oscillations 
gradually decay to a uniform bulk proflle within 5 — 6 a away from the wall (not shown). 
The contact density decreases slightly for larger values of the applied force. A correlation 
between the fluid structure near the wall and the slip length will be examined in the next 
section. 

The solution of the Navier-Stokes equation for the force driven flow with slip boundary 
conditions at the conflning parallel walls, v(±/i) = Vg-, is given by [l] 



v(^) = ^ (/^' - z^) + Vg 
2,(1 



(6) 
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FIG. 2: (Color online) Averaged velocity profiles, (v)r/(T, for the indicated values of the applied 
force per fluid molecule for the rigid walls (left) and the thermal walls with K = 600e/(T^ (right). 
The wall-fluid interaction energy is fixed to e„f/e = 0.9. The solid lines represent a parabolic fit to 
the data. The dashed lines indicate the location of liquid/solid interfaces at z = ± 11.79(7. Vertical 
axes denote the position of the fee lattice planes at z = ±12.29 o". 

where 2h is the distance between the walls and /i is the shear rate independent viscosity. The 
slip length is defined as an extrapolated distance relative to the position of the liquid/solid 
interface where the tangential component of the fluid velocity vanishes 



(7) 



Figure[2] shows representative velocity profiles in steady-state flow for three different values 
of the external force fx and fixed wall-fluid interaction energy £:„f/£: = 0.9. The data for 
thermal and rigid walls are presented only in half of the channel because of the symmetry 
with respect to the mid-plane of the fluid phase, i.e. v{z) = v{—z). Fluid velocity profiles 
are well fitted by a parabola with a shift by the value of the slip velocity, as expected from 
the continuum predictions [see Eq. (^]. The simulation results presented in Fig. [2] show that 
slip velocity Vs increases with the applied force. The degree of slip depends on the wall 
stiffness and the interfacial shear rate. The fluid slip velocity is larger for the thermal walls 
and small forces fx ^ 0.012 e/a. By contrast, rigid walls produce more slippage for the large 
value of the external force fx = 0.024 er/a. In both cases, the slip velocity is greater than the 
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le channel and near the walls. The 



261], ensuring laminar flow conditions 



difference between the fluid velocities at the center o: 
upper bound for the Reynolds number is Re «i 10 
throughout. 



B. Effect of thermal wall roughness on slip length 

Slip boundary conditions for a fluid flow past atomically smooth rigid walls are determined 
by the molecular-scale surface roughness due to the wall atoms flxed at their equilibrium 
lattice sites. The thermal fluctuations of the wall atoms, being unavoidable in real surfaces, 
modify the effective coupling between liquid and solid phases. Depending on the stiffness 
of the surface, thermal walls can either reduce slip (due to the deep penetration of the wall 
atoms into the fluid phase) or increase slip because of the reduction of the surface induced 
structure in the adjacent fluid layer. In this section, the shear rate dependence of the slip 
length is investigated in a wide range of values of the parameter k that satisfy the Lindemann 
criterion for melting (see Table [T]). 



In the previous MD study on shear flow near solids by Thompson and Robbins [I4], it 
was shown that the thermal surface roughness reduces the degree of ordering in the adjacent 
fluid layer. The thermal fluctuations of the wall atoms produced slip lengths of about 0.5 a 
larger than their values for the rigid walls. In a wide range of wall-fluid interaction energies 
0.2 ^ £:„f/£: ^ 25, the slip length was found to be rate-independent and less than 3.5 a. In 
our simulations, the surface potential is less corrugated because of the higher wall density; 
and, therefore, the effect of thermal surface roughness on the slip length is greater. 

The variation of the slip length with increasing shear rate for different values of the spring 
stiffness n and e^f/e = 0.9 is presented in Fig.[3l The data for smooth rigid walls, fltted by a 
straight line, are also shown in Fig. [3] for comparison with the results for the thermal walls. 
In the range of accessible shear rates, the slip length is larger for stiffer thermal walls. The 
surface becomes effectively smoother because the average displacement of the wall atoms 
with respect to their equilibrium sites is reduced at larger values of k. A similar reduction 
in slip velocity for the soft thermal wall atoms was reported in recent MD simulations of 



thin fllms of hexadecane 211]. For soft walls with K, = A00e/a^, the wall atoms penetrate 
deeper into the fluid phase, which makes the slip length smaller and weakly dependent on 
the shear rate. For k = 600 e/o"^, the slip length increases slightly at low shear rates and 
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FIG. 3: (Color online) Slip length, Lg/a, as a function of the shear rate at the interface for 
£^f/e = 0.9. The values of the spring constant k for the thermal walls are listed in the inset. The 
dashed lines represent the best fit to the data for the thermal walls with K = 1600e/(7^ (v) and 
the rigid fee walls (o). The slip length for rough rigid walls with 0.07 cj ([>). The solid 

curves are a guide for the eye. 

then saturates at 'jr > 0.063, where it becomes smaller than Ls for the rigid walls. This 
behavior is consistent with the results for the fluid velocity profiles presented in Fig. [2] for 
the thermal and rigid walls. 

In the case of the largest spring constant k = 1600 e/cr^, the slip length increases mono- 
tonically with the shear rate and its dependence can also be fitted well by a straight line. 
The slope of the fitted line is slightly larger than one for the rigid walls (see Fig.[3]). For 
a finite spring stiffness, a small downward curvature appears at 7r > 0.05 because of the 
higher temperature and, as a consequence, larger mean displacement of the wall atoms. The 
maximum increase in temperature of the wall atoms and the adjacent fluid layer is about 
10% at the highest shear rates reported in Fig.O The fluid temperature at the center of 
the channel increases up to 1.106 e/fc^ at the highest 7. The upper bound for the shear 
rate is determined by the maximum shear stress the liquid/solid interface can support {26 ]. 
The simulation results for the thermal walls presented in Fig. [3] demonstrate that the spring 
stiffness in the model of harmonic oscillators can be an important factor in determining the 
degree of slip in the rate-dependent regime. 
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FIG. 4: (Color online) (a) Behavior of the slip length Lg/cr as a function of the inverse value of the 
in-plane structure factor, l/S'(Gi), evaluated at the first reciprocal lattice vector Gi = (9.04(T^^, 0). 
The wall-fluid interaction energy is £^ile = 0.9. (b) The same data for the indicated values of the 
spring constant are replotted versus [5(Gi) pcO^\'^ . 

It is interesting to note that the surface roughness due to immobile wall atoms with 
random displacement of only a fraction of a molecular diameter significantly reduces slip 
length and leads to a slight upward curvature in the rate dependence (see Fig. [3]). The rough 
surfaces were constructed by fixing the instantaneous positions of initially equilibrated wall 
atoms with the spring stiffness k = 600 e/a^ in the absence of the flow. The parabolic velocity 
profiles for different values of the applied force were computed for the same realization of 
disorder. These random perturbations of the surface potential lead to the difference in the 
slip length of about 7o" in comparison with its values for the thermal walls with the spring 
stiffness k= 600 e/a^. 

The degree of slip at the liquid/solid interface correlates well with the amount of the 



surface induced order in the adjacent fluid layer [1^, |20|, l23|, |36|. The in-plane fluid struc- 
ture factor is deflned as 5'(k) = 1/Ne \ J2j e^^'^^\^, where = {xj, yj) is the two-dimensional 
position vector of the j-th molecule and the sum is taken over Ni molecules within the flrst 
layer. The effect of periodic surface potential on the structure of the adjacent fluid becomes 
more pronounced at the reciprocal lattice vectors. In the previous MD study j26| for sim- 
ilar parameters of the wall and fluid phases, it was shown that the slip length scales as 
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Ls~(Ti/S'(Gi) Pc)") where Gi is the the first reciprocal lattice vector in the flow direction, 
Ti is temperature of the first fluid layer and a = 1.44 ±0.10. This scaling relation was found 
to hold in a wide range of shear rates and wall-fluid interactions for atomically smooth rigid 
walls and incommensurate structures of the liquid/solid interface {26]. 

The correlation between the inverse value of the fluid structure factor evaluated at the 
first reciprocal lattice vector Gi = (9.04 o""^, 0) and the slip length is presented in Fig.|l](a). 
Except for the rough rigid walls, the surface induced structure in the first fluid layer 5'(Gi) is 
reduced at higher shear rates and smaller values of k. The slip length increases approximately 
linearly with l/S'(Gi) for the rigid and stiff walls with k, ^ 1200 e/a^. A gradual transition 
to a weak dependence of the slip length on S'(Gi) is observed upon reducing the wall stiffness. 
In Figure 111(b) the same data for the slip length are replotted as a function of the inverse 
product [5'(Gi) pccr^]~^, where pc is a contact density of the first fluid layer. Although the 
contact density decreases slightly with increasing the slip velocity, the functional form of 
the slip length is similar in both cases [see Figs.|4](a) -|4](b)]. The results shown in Fig.H] 
indicate that the dependence of the slip length on the fluid structure of the first layer is 
less pronounced for k, ^ 800 e/a"^ due the thermal surface roughness in the rate-dependent 



26| holds in the presence of the 



regime. Whether the scaling relation for the slip length 
thermal roughness for different wall-fluid interaction energies will be the subject of the future 
research. 



IV. RESULTS FOR PERIODICALLY CORRUGATED WALLS 

Next, the results for the rate dependence of the slip length are compared for atomically 
smooth rigid walls and periodically roughened surfaces. The periodic surface roughness of 
the upper and the lower walls was modeled by introducing a vertical offset to the positions 
of the wall atoms Az{x) = a sin(27rx/A) with the wavelength A = 4.17 a. In this part of the 
study, the wall atoms are rigidly fixed with respect to their equilibrium sites. To properly 
compare the results for atomically smooth and roughened surfaces, both the local shear rate 
and the slip length were estimated from a parabolic fit of the velocity profiles at the same 
location of the interface, i.e. z = ± 11.79 a. A weaker wall-fiuid interaction, e^ile = 0.5, was 
chosen to obtain larger values of in the absence of the imposed corrugation, since it is 
expected that the surface roughness strongly reduces the slip length 
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FIG. 5: (Color online) Variation of the slip length Lg/cr as a function of the local shear rate 
for Poiseuille flows and £^{/e = 0.5. The wavelength and amplitudes of the wall modulation are 
tabulated in the insets. The local shear rate and the slip length are extracted from a parabolic fit 
of the velocity profiles at z = ± 11.79 cr. 

The dynamic response of the slip length with increasing shear rate is presented in Fig. [5] 
for atomically smooth and periodically corrugated walls. The data for the flat walls (a = 0) 
are the same as in Ref. j^]- At low shear rates, the slip length L° (defined by the leftmost 
point on each curve in Fig. [5]) decreases monotonically with increasing the amplitude a. 
The no-slip boundary condition is achieved for a > 0.3 a. This behavior for a similar 
geometry and interaction parameters was examined in detail in the previous paper 28|. 
A direct comparison between continuum analysis and MD simulations showed that there 
is an excellent agreement between the velocity profiles and the slip length for the large 
wavelengths A > 20cr and small values of a/A < 0.05. The continuum results overestimate 
the slip length when A approaches molecular dimensions 



At higher shear rates, the slope of the rate-dependent slip length is gradually reduced with 
increasing the amplitude of the surface corrugation (see Fig.E]). For the largest amplitude 
a = 0.3 a, the slip length weakly depends on shear rate and its magnitude becomes smaller 
than the molecular diameter. As apparent from the set of curves shown in Fig.[5l the same 
value of the slip length can be obtained by increasing simultaneously the amplitude of the 
surface corrugation and the shear rate. Analogous behavior of the slip length was observed 
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experimentally for flows of Newtonian liquids past surfaces with variable nano-roughness 31 1 ■ 
We note, however, that the MD simulations of simple fluids reported in this study do not 
show any threshold in the rate dependence of the slip length for the amplitudes of the surface 
corrugation a ^ 0.3 a. 



V. SUMMARY 



In this paper the effect of molecular-scale surface roughness on the slip length in a flow 
of simple fluids was studied by molecular dynamics simulations. The parabolic flt of the 
steady-state velocity proflles induced by a constant force was used to deflne the values of 
interfacial shear rate and slip length. For atomically smooth rigid surfaces and weak wall- 
fluid interactions, the slip length increases approximately linearly with the shear rate. Three 
types of surface roughness were considered: thermal, random and periodic. 

The thermal surface roughness due to flnite spring stiffness of the wall atoms signiflcantly 
modifies the slip behavior. The large penetration of the wall atoms into the fluid phase ob- 
served for soft walls causes weak rate dependence of the slip length below its values for 
atomically smooth rigid walls. Increasing the wall stiffness produces effectively smoother 
surfaces and leads to the linear rate dependence of the slip length. Periodically and ran- 
domly corrugated rigid surfaces, with the amplitude below the molecular diameter, strongly 
reduce the slip length and its shear rate dependence. These findings open perspectives for 
modeling complex systems with combined effects of surface roughness, wettability and rate 
dependency. 
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